library(DESeq2)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pheatmap)
library(outliers)
library(dplyr)
library(plotly)
library(grid)
library(futile.logger)
library(VennDiagram)
options(width=120)

В ходе данной домашней работы понадобятся те же файлы, что и в лекции: “GSE89225_illumina_counts.csv”, “conditions.csv”, “human_mart.txt”. Для начала убедимся в том, что мы можем эти файлы прочитать. И посмотрим, что в них находится.

Неуместного аутлайера с условным обозначением treg_NBP_patient3 я вычислила, но немного позднее по ходу этого файла. Но, хотя Читатель узнает его имя позднее, я, раз знаю немного больше, удалю его из данных прямо сейчас сразу после считывания файла.

counts <- read.csv("GSE89225_Illumina_counts.csv", row.names=1)
conditions <- read.csv("conditions.csv", row.names=1)
mart <- read.table("human_mart.txt", sep="\t", header=1, check.names = F)

# а вот и его ненастоящее имя treg_NBP_patient3

counts <- select(counts, -treg_NBP_patient3)

conditions <- conditions[-12,]

# а это несколько строчек, чтобы взглянуть на струкутру датафреймов 
print(counts[1:6, 1:2])
##                 tconv_tumor_patient10 cd177negtreg_tumor_patient10
## ENSG00000000003                   510                          935
## ENSG00000000005                    25                            1
## ENSG00000000419                   613                          818
## ENSG00000000457                   269                          587
## ENSG00000000460                    55                          142
## ENSG00000000938                    36                           69
dim(counts)
## [1] 63677    33
head(conditions)
##                                            tissue                               cells
## tconv_tumor_patient10        tissue: breast tumor cell type: Conventional CD4 T cells
## cd177negtreg_tumor_patient10 tissue: breast tumor       cell type: Regulatory T cells
## cd177postreg_tumor_patient10 tissue: breast tumor       cell type: Regulatory T cells
## tconv_tumor_patient1         tissue: breast tumor cell type: Conventional CD4 T cells
## treg_tumor_patient1          tissue: breast tumor       cell type: Regulatory T cells
## tconv_NBP_patient2                    tissue: NBP cell type: Conventional CD4 T cells
dim(conditions)
## [1] 33  2
head(mart)
##           Gene ID Associated Gene Name Gene type EntrezGene ID
## 1 ENSG00000252303            RNU6-280P     snRNA            NA
## 2 ENSG00000281771                Y_RNA  misc_RNA            NA
## 3 ENSG00000281256         RP11-222G7.2   lincRNA            NA
## 4 ENSG00000283272      Clostridiales-1      sRNA            NA
## 5 ENSG00000280864        RP11-654C22.2   lincRNA            NA
## 6 ENSG00000280792        RP11-315F22.1   lincRNA            NA
dim(mart)
## [1] 64101     4

RNA-seq

Немного слов учителя об РНК-секвенировании:

Rna-seq steps:

  • Изоляция РНК
  • Rna selection / depletion
  • Фрагментация
  • Синтез кДНК
  • Секвенирование

Rna selection / depletion:

  • вся РНК
  • тянем за поли(А)-хвосты (только мРНК)
  • удаляем рибосомальную РНК (ribo-zero kit)
  • таргетное секвенирование

Why Rna-seq?

  • Не ограничены существующей сборкой и дизайном микрочипа
  • Низкий фоновый сигнал
  • Точность позволяет смотреть экспрессию отдельных изоформ

Sanity checks

###…и необходимой чистоте данных

Нужно всегда проверять длины библиотек и количество rRNA reads, которые оказались в библиотеке. Количество ридов можно проверять после выравнивания или после квантификации.

# Повторим сказанное на языке машины

proteinCoding <- mart[mart[, 3] == "protein_coding", ]
rRNA <- mart[mart[, 3] == "rRNA", ]

pcCounts <- counts[rownames(counts) %in% as.character(proteinCoding[, 1]), ]
rrnaCounts <- counts[rownames(counts) %in% as.character(rRNA[, 1]), ]

sampleCount <- ncol(counts)
toPlot <- data.frame(
  sample=rep(colnames(counts), 3),
  value=c(colSums(counts) - colSums(pcCounts) - colSums(rrnaCounts), 
          colSums(pcCounts), 
          colSums(rrnaCounts)),
  type=c(rep("other", sampleCount), 
         rep("protein coding", sampleCount),
         rep("rrna", sampleCount))
)

plot <- ggplot(data=toPlot, aes(x=sample, y=value, fill=type)) +
  geom_bar(stat="identity") + theme_bw() + 
  theme(axis.text.x = element_text(angle=90, vjust=0.5))
plot

DESeq2

DESeq2 и правда хорош и, возможно, даже красив. Тут и дифференциальная экспрессия, и нормализации, и PCA-plots.

# исправилась, была неправа: вот теперь два DESeq-объекта

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = conditions,
                              design = ~ tissue + cells)

dds <- dds[rowSums(counts(dds)) > 20, ]
dds <- DESeq(dds)
vst_dds <- vst(dds)
counts.norm <- assay(vst_dds)

# а вот и тот самый вышеупомянутый outlier, встречаем: 

outlier(counts.norm)
##        tconv_tumor_patient10 cd177negtreg_tumor_patient10 cd177postreg_tumor_patient10         tconv_tumor_patient1 
##                     18.31254                     19.06682                     18.75782                     19.50822 
##          treg_tumor_patient1           tconv_NBP_patient2            treg_NBP_patient2         tconv_tumor_patient2 
##                     19.42026                     18.74975                     18.97194                     18.68064 
##  cd177negtreg_tumor_patient2  cd177postreg_tumor_patient2           tconv_NBP_patient3         tconv_tumor_patient3 
##                     18.86007                     19.00433                     18.76385                     19.09414 
##  cd177negtreg_tumor_patient3  cd177postreg_tumor_patient3           tconv_NBP_patient4            treg_NBP_patient4 
##                     18.73277                     19.78128                     18.74858                     19.24072 
##         tconv_tumor_patient4  cd177negtreg_tumor_patient4  cd177postreg_tumor_patient4          treg_tumor_patient5 
##                     18.51430                     19.73984                     19.74446                     18.87674 
##         tconv_tumor_patient5         tconv_tumor_patient6  cd177negtreg_tumor_patient6  cd177postreg_tumor_patient6 
##                     18.84609                     18.85645                     18.68582                     18.93730 
##           tconv_NBP_patient7            treg_NBP_patient7           tconv_NBP_patient8            treg_NBP_patient8 
##                     18.55074                     18.03944                     18.73227                     18.73946 
##          treg_tumor_patient8         tconv_tumor_patient8           tconv_NBP_patient9            treg_NBP_patient9 
##                     18.43248                     18.40272                     18.90248                     18.56628 
##         tconv_tumor_patient9 
##                     18.94650
max(outlier(counts.norm)) # а это treg_NBP_patient3
## [1] 19.78128
# Вот такой еще

dds_cells <- DESeqDataSetFromMatrix(countData = counts,
                              colData = conditions,
                              design = ~ cells + tissue)

dds_cells <- dds_cells[rowSums(counts(dds_cells)) > 20, ]
dds_cells <- DESeq(dds_cells)
vst_dds_cells <- vst(dds_cells)
counts.norm_cells <- assay(vst_dds_cells)

Поглядим на PCA без выброса?

pca_data_tissue <- prcomp(t(counts.norm))
percents_tissue <- pca_data_tissue$sdev^2 / sum(pca_data_tissue$sdev^2)
to_plot_tissue <- t(counts.norm) %*% pca_data_tissue$rotation

pca_data_cells <- prcomp(t(counts.norm_cells))
percents_cells <- pca_data_cells$sdev^2 / sum(pca_data_cells$sdev^2)
to_plot_cells <- t(counts.norm_cells) %*% pca_data_cells$rotation

gdata_tissue <- data.frame(
  x=to_plot_tissue[, 1],
  y=to_plot_tissue[, 2],
  tissue=conditions[, 1],
  name=rownames(conditions)
)

gdata_cells <- data.frame(
  x=to_plot_cells[, 1],
  y=to_plot_cells[, 2],
  cells=conditions[, 2],
  name=rownames(conditions)
)

# немного ненавязчивого интерактива от plotly

ggplotly(ggplot(data=gdata_tissue, aes(x=x, y=y, color=tissue, shape=tissue, text=name)) +
  geom_point(size=3) + theme_bw()  +
  xlab(paste0("PC", 1, ": ", formatC(100 * percents_tissue[1], digits=4), "%")) +
  ylab(paste0("PC", 2, ": ", formatC(100 * percents_tissue[2], digits=4), "%")))
ggplotly(ggplot(data=gdata_cells, aes(x=x, y=y, color=cells, shape=cells, text=name)) +
  geom_point(size=3) + theme_bw()  +
  xlab(paste0("PC", 1, ": ", formatC(100 * percents_cells[1], digits=4), "%")) +
  ylab(paste0("PC", 2, ": ", formatC(100 * percents_cells[2], digits=4), "%")))
ggplotly(plotPCA(vst_dds, intgroup=c("tissue", "cells"), ntop = 500) + theme_bw())

Differential expression

“Давайте посмотрим, как выглядят результаты дифференциальной экспрессии и отсортируем их по статистике”.

res <- results(dds)
res <- res[order(res[,4]),]

res_cells <- results(dds_cells)
res_cells <- res_cells[order(res_cells[, 4]), ]

Наконец исправила Volcano. А знаешь, что было не так? В объекте dds_cells на 135 строчке у меня была русская с!!!!!!!!!!!!!!!!!!!!

genes_tissue <- rownames(res)
genes_cells <- rownames(res_cells)

gdata_tissue <- data.frame(
  log_fold_change=res$log2FoldChange,
  p_adj=res$padj,
  n = "Treg vs Tconv"
)

gdata_tissue <- na.omit(gdata_tissue)

gdata_cells <- data.frame(
  log_fold_change=res_cells$log2FoldChange,
  p_adj=res_cells$padj,
  n = "Breast tumor vs Normal breast tissue"
)

gdata_cells <- na.omit(gdata_cells)

# я все это делаю еще и потому, что потом мне отдельно понадобятся датафреймы для клеток и тканей

# но для volcano, хотя я и планировала сначала использовать grid_arrange (почему и создавала два gdata для клеток и для тканей), в конечном итоге я использую слитый воедино датафрейм general_gdata. Потому как в последующем применю facet_grid.

general_gdata <- rbind(gdata_tissue, gdata_cells)
general_gdata <- na.omit(general_gdata)

# создадим порог значимости

general_gdata$threshold <- general_gdata$p_adj <= 0.01
general_gdata$significant <- sapply(general_gdata$threshold, function(x) ifelse(x == TRUE, "Significant", "Not Significant"))
general_gdata <- na.omit(general_gdata)

# Вы только посмотрите, какая красота!

ggplot(general_gdata,aes(x=log_fold_change, y=-log10(p_adj), col = significant)) +
       facet_grid(. ~ n) +
       geom_point(size=1) +
       theme_bw() +
       scale_colour_manual(values = c("black", "red")) + 
       labs(x = "Log fold change", y ="Adjusted p.value")+
       geom_hline(yintercept=2, linetype = 2, size = 1, colour = 'red')

Поскольку я задержалась с домашкой, я решила, а может, еще что-нибудь почитать и как-нибудь поиграть с построением volcano? Прием, которым я воспользовалась ниже, направлен на более тонкое выявление дифференциальной экспрессии, основанное и на значении log2 fold change между двумя уровнями факторной переменной (уровнем экспрессии генов в разных клетках и в разных тканях) и на значении p-adj (p-val после поправки BH (Benjamini & Hochberg (1995)). Из полученного графика мы видим, что хотя различия в экспресии по показателю log2 fold change между некоторыми генами и выявляются, они остаются тем не менее незначимыми (черным цветом).

general_res <- rbind(res, res_cells)

p_adj <- c(general_res$padj)

log_fold_change <- c(general_res$log2FoldChange)

this_frame <- general_gdata[,-c(4,5)]

this_frame_green <- subset(this_frame, log_fold_change < -1 & p_adj < .01) # определим зеленую зону

this_frame_green <- cbind(this_frame_green, rep(1, nrow(this_frame_green)))

colnames(this_frame_green)[4] <- "Color"

this_frame_black <- subset(this_frame, (log_fold_change >= -1 & log_fold_change <= 1) | p_adj >= .01) # определим плохую черную зону

this_frame_black <- cbind(this_frame_black, rep(2, nrow(this_frame_black)))

colnames(this_frame_black)[4] <- "Color"

this_frame_red <- subset(this_frame, log_fold_change > 1 & p_adj < .01) # определим красную зону

this_frame_red <- cbind(this_frame_red, rep(3, nrow(this_frame_red)))

colnames(this_frame_red)[4] <- "Color"

frame_total <- rbind(this_frame_green, this_frame_black, this_frame_red)

frame_total$Color <- as.factor(frame_total$Color)

genes <- rownames(general_res)

##Сконструируем объект volcano

my_another_volcano_hello <- ggplot(data = frame_total, aes(x = log_fold_change, y = -log10(p_adj),col= Color)) +
  geom_point(alpha = 0.5, size = 1) + theme_bw() + theme(legend.position = "none") +
  geom_hline(yintercept = 2, colour = "black", linetype = 2, size = 1) + 
  scale_color_manual(values = c("green", "black", "red")) +
  facet_grid(. ~ n) +
  labs(y="Adjusted p.value", x="Log fold change")

ggplotly(my_another_volcano_hello)

А теперь взглянем на горячую карту. Горячая карта покажет горячие точки. Вдруг эти точки о чем-то говорят. (Шучу, - конечно же, они говорят)

# процедуры извлечения генов из пасвея

kkeys <- keys(org.Hs.eg.db, keytype="ENSEMBL")
goAnno <- AnnotationDbi::select(org.Hs.eg.db, keys=kkeys, 
                                keytype="ENSEMBL", columns=c("GOALL", "ONTOLOGYALL", "SYMBOL"))
goAnno <- tbl_df(goAnno)
goAnno <- filter(goAnno, GOALL=="GO:0007159")
# or you can pick ENTREZ, or SYMBOL, or whatever you want
genesToVisualise <- goAnno$ENSEMBL

# несколько шагов перед тем, как предложить что-то изящной функции

to_visualise_zzzz <- counts.norm[rownames(res), order(conditions[, 2])]

subsetting <- rownames(subset(res, rownames(res) %in% genesToVisualise))

nearly_want_to_see <- to_visualise_zzzz[subsetting,]

exactly_want_to_see <- t(apply(nearly_want_to_see, 1, function(r) {(r - min(r)) / (max(r) - min(r))}))

# О, как она изящна

pheatmap(exactly_want_to_see, 
         show_rownames = F, cluster_rows = F,
         cluster_cols=F,
         annotation_col = conditions,
         main = "GO:0007159: leukocyte cell-cell adhesion")

Построим же наконец диаграмму Венна

for_set_1 <- rownames(res[which(res$padj < 0.01),])
for_set_2 <- rownames(res_cells[which(res_cells$padj < 0.01),])

wanted_sected <- intersect(for_set_1, for_set_2)

length(wanted_sected) # 84
## [1] 84
venn_plot <- draw.pairwise.venn(length(for_set_1), length(for_set_2), length(wanted_sected), category = c("Tregs vs Tconv", "Tissues"), fill = c("blue", "pink"), ext.line.lty = "dashed")